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Abstract 

We describe a new, free and open source semi-analytic model of galaxy formation, Galacticus. The Galacticus 
model was designed to be highly modular to facilitate expansion and the exploration of alternative descriptions of 
key physical ingredients. We detail the Galacticus engine for evolving galaxies through a merging hierarchy of dark 
matter halos and give details of the specific implementations of physics currently available in Galacticus. Finally, we 
show results from an example model that is in reasonably good agreement with several observational datasets. We 
use this model to explore numerical convergence and to demonstrate the types of information which can be extracted 
from Galacticus. 
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1. Introduction 

The physics of galaxy formation is rich and complex, 
and has provided a challenge for theorists to understand 
and explain ever since it became clear that the Uni- 



verse is filled with galaxies (Shapley and Curtis 1921 1. 



Various theoretical tools have been brought to bear on 
this problem, ranging from direct analytic reasoning to 
large scale numerical N-body simulations. An inter- 
mediate and highly successful approach is that known 
as "semi-analytic" galaxy formation modeling. In this 
approach the numerous complex non-linear physics in- 
volved are solved using a combination of analytic ap- 
proximations and empirical calibrations from more de- 
tailed, numerical solutions. Historically, such models 
were first contemplated by | White and Rees ( 1978[> and 



have since been developed further, notably by 



White 



|and Frenk| ( fl99"TT >, [Kauffmann et al.| ( fl993l >, [Cole et al. 
pOOO) , |Hatton et aL[ ( pOOl} , |Monaco et aE[ ( [2007] ), 
|Somerville et al. ( 2008[ l among others. Models of this 
type aim to begin with the initial state of the Universe 
(specified shortly after the Big Bang) and apply physi- 
cal principles to determine the properties of galaxies in 
the Universe at later times, including the present day. 
Typical properties computed include the mass of stars 
and gas in each galaxy, broad structural properties (e.g. 
radii, rotation speeds, geometrical shape etc.), dark mat- 
ter and black hole contents, and observable quantities 
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such as luminosities, chemical composition etc. As 
with N-body/hydrodynamical simulations, the degree of 
approximation varies considerably with the complexity 
of the physics being treated, ranging from precision- 
calibrated estimates of dark matter merger rates to em- 
pirically motivated scaling functions with large param- 
eter uncertainty (e.g. in the case of star formation and 
feedback). 

The primary advantage of the semi-analytic approach 
is that it is computationally inexpensive compared to N- 
body/hydrodynamical simulations. This facilitates the 
construction of samples of galaxies orders of magnitude 
larger than possible with N-body techniques and for the 
rapid exploration of parameter space (Henriques et al., 
[20091 iBeiisbn and Bower) |2010b| [Bower et al.| |20ioj l 
and model space (i.e. adding in new physics and as- 
sessing the effects). The primary disadvantage is that 
they involve a greater degree of approximation. The ex- 
tent to which this actually matters has not yet been well 
assessed. Comparison studies of semi-analytic vs. N- 
body/hydro calculations have shown overall quite good 
agreement (at least on mass scales well above the reso- 
lution limit of the simulation) but have been limited to 
either simplified physics (e.g. hydrodynamics and cool- 
ing only;|Benson et al.|2"00T[|Yoshida et al.|2002[|Heliy" 
|et al.|2003[|Lu et al.|2010) > or to simulations of individ- 
ual galaxies ( Stringer et al.||2010[ l. 

In this work, we describe a new, free and open 
source semi-analytic model, Galacticus. The Galacti- 
cus model solves the physics describing how galaxies 
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evolve in a merging hierarchy of dark matter halos in a 
cold dark matter universe. Galacticus has much in com- 
mon with other semi-analytic models, such as the range 
of physical processes included and the type of quanti- 
ties that it can predict, but has some key distinguishing 
features. In designing Galacticus our main goal was to 
make the code flexible, modular and easily extensible. 
Much greater priority was placed on making the code 
easy to use and to modify than on making it fast. We 
believe that a modular and extensible nature is crucial 
as galaxy formation is an evolving science. In particu- 
lar, key design features are: 

Extensible methods for all functions: Essentially all 
functions within Galacticus are designed to be ex- 
tensible, meaning that anyone can write their own 
version and easily insert it into Galacticus. For 
example, suppose an improved functional form for 
the cold dark matter (CDM) halo mass function is 
derived. A user of Galacticus can simply write a 
short module conforming to a specified template 
that computes this mass function and which in- 
cludes a short directive in the code which explains 
to Galacticus's build system how to incorporate 
this module. A recompile of the code will then in- 
corporate the new function. The user is absolved 
from having to understand the details of the in- 
ner workings of the code, instead simply being re- 
quired to conform to a standard interface. 

Extensible components for tree nodes: The basic 
structure in Galacticus is a merger tree, which 
consists of a set of linked tree nodes which 
have various properties. Galacticus works by 
evolving the nodes forward in time subject to a 
collection of differential equations and other rules. 
Each node can contain an arbitrary number of 
components. For example, a component may be a 
dark matter halo, a galactic disk, a black hole etc. 
Each component may have an arbitrary number 
of properties (some of which may be evolving, 
others of which can be fixed). Galacticus makes it 
easy to add additional components. For example, 
suppose that a user wanted to add a "stellar halo" 
component (consisting of stars stripped from 
satellite galaxies). To do this, they would write 
a module which specifies the following for this 
component: 

• Number of properties; 

• Interfaces to set and get property values and 
rates of change; 



• "Pipes" which allow for flows of 
mass/energy/etc. from one component 
to another (and which facilitate interaction 
between components without the need for 
knowledge of the specific implementation of 
any component); 

• Functions describing the differential equa- 
tions which govern the evolution of the prop- 
erties; 

• Functions describing how the component re- 
sponds to various events (e.g. the node be- 
coming a satellite, a galaxy-galaxy merger, 
etc.); 

• Auxiliary routines for handling outputs etc. 

Short directives embedded in this module explain 
to the Galacticus build system how to incorporate 
the new component. A recompile will then build 
the new component into Galacticus. Typically, a 
new component can be created quickly by copy- 
ing an existing one and modifying it as necessary. 
Furthermore, multiple implementations of a com- 
ponent are allowed. For example, Galacticus con- 
tains a component which is a Hernquist spheroid. 
One could add a de Vaucouler's spheroid compo- 
nent and an input parameter then allows one to sim- 
ply select which implementation will be used in a 
given run. 

Centralized ODE solver: Galacticus evolves nodes 
in merger trees by calling an ODE solver which in- 
tegrates forward in time to solve for the evolution 
of the properties of each component in a node. This 
means that it is not necessary to provide explicit 
solutions for ODEs (in many cases such solutions 
are not available anyway) and time-stepping is au- 
tomatically handled to achieve a specified level of 
precision. The ODE solver allows for the evolu- 
tion to be interrupted. A component may trigger 
an interrupt at any time and may do so for a num- 
ber of reasons. A typical use is to actually cre- 
ate a component within a given node — for example 
when gas first begins to cool and inflow in a node 
a galactic disk component may be created. Other 
uses include interrupting evolution when a merg- 
ing event occurs. 

To summarize, the Galacticus code is therefore 
highly modular. Every part of it consists of a simple 
and well-defined interface into which an alternative im- 
plementation of a calculation can easily be added. Sim- 
ilarly, the physical description of galaxies is extremely 
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flexible. Each galaxy has a set of components which 
can be created/destroyed as needed, each of which has 
a set of properties. These components are also modu- 
lar, and any such module can be trivially replaced by 
another that performs the calculations differently if re- 
quired. The actual formation and evolution of galaxies 
is treated by simply defining a set of differential equa- 
tions for each galaxy. Then, these are all fed in to the 
ODE solver which evolves them to a specified accuracy. 

The goal of this article is to describe the technical 
and physical implementation of Galacticus. De- 
tailed examination of its scientific predictions and 
their consequences will be deferred to a future work. 
The remainder of this article is arranged as follows. 
In ^2] we describe the technical implementation of 
Galacticus, while in ^3] we describe the current spe- 
cific implementation of galactic components, and in 
^4] the current specific implementation of numerous 
physical processes and properties. In ^5] we show the 
results of an example calculation. Finally, in ^6] we 
summarize the information presented in this article. 
The Galacticus model is available for download from 
http : //sites . google . com/site/galacticusmodel 
A full manual documenting both the physics and tech- 
nical implementation of Galacticus can also be found 
on the same website. 

2. Technical Implementation 

In the following section we describe the technical and 
practical aspects of Galacticus. The Galacticus code is 
free and open source and was designed to depend only 
on free and open source compilers and libraries to max- 
imize its portability. In addition to the GNU compilers 
Galacticus requires the GNU Scientific Library (GSL) 
and the FGSL wrapper for GSL, the FoX XML parser 
and the HDF5 library. Additionally, the analysis codes 
provided with Galacticus make extensive use of Perl 
and PDL although Galacticus output can be just as eas- 
ily analyzed with other tools. 

2.1. Running 

The behavior of Galacticus is controlled by a set of 
parameters which are given in a file specified as a com- 
mand line argument. (Galacticus will in fact provide 
sensible default values for all parameters if no param- 
eter file is specified, or if some parameters are miss- 
ing from the file). The parameter file is an XML file, 
which allows it to be generated easily using a variety of 
pre-existing XML tools and libraries. Scripts are pro- 
vided with Galacticus to generate parameter files and 



run grids of models which span a range of parameter 
values. 

When run, Galacticus proceeds by performing a set 
of predefined tasks in order until all tasks are done. As 
with all other aspects of Galacticus, these tasks are 
modular and extensible, allowing new tasks to be added 
in as desired. Typically, however, the tasks will consist 
of some initialization, followed by creation and evolu- 
tion of one or more merger trees. 

2.2. Output 

Galacticus stores its output in an HDF5 file which 
allows this output to be viewed and manipulated using a 
variety of tools in addition to the standard HDF5 C API 
including: 

HDFView A graphical viewer for exploring the con- 
tents of HDF5 files; 

HDF5 Command Line Tools A set of tools which can 
be used to extract data from HDF5 files (h5dump 
and h51s are particularly useful); 

C++ and Fortran 90 APIs Allow access to and ma- 
nipulation of data in HDF5 files from within C++ 
and Fortran90 codes; 

h5py A Python interface to HDF5 files; 

IPDL::IO::HDF5l A Perl interface to HDF5 files. 

2.2.1. Output Datasets 

In addition to the properties of galaxies and their 
associated dark matter halos in all merger trees, each 
Galacticus output file stores a full record of all param- 
eter values (either input or default) that were used for 
the particular run. These can easily be extracted to an 
XML file suitable for re-input into Galacticus. The 
output also contains a record of the Galacticus version 
used for this model, storing the major and minor version 
numbers, and the revision number along with the time 
at which the model was run. 

Optionally, Galacticus will compute and store vol- 
ume averaged properties of the entire galaxy popula- 
tion at a set of snapshot times. Currently, the proper- 
ties stored are star formation rate density, stellar mass 
density, interstellar medium (ISM) gas density and the 
density in resolved dark matter halos, all as a function 
of cosmic time. 

Galaxy data can be output at one or more snapshot 
times, specified as input parameters to Galacticus. At 
each output time, each merger tree is stored separately 
within the HDF5 file to facilitate easy access. Each such 
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merger tree group contains all data on a single merger 
tree, and consists of a collection of datasets each of 
which lists a property of all nodes in the tree which ex- 
ist at the output time. Additionally, a weight (in Mpc~ 3 ) 
which should be assigned to this tree (and all nodes in 
it) to create a volume-averaged sample is stored. 

Optionally, a mass accretion history (i.e. mass as a 
function of time) for the main branctrl in each merger 
tree can be output, which lists the mass of the main 
branch as a function of time. 

Finally, Galacticus can output the full structure of 
merger trees prior to any evolution. This is useful to 
permit the same trees to be used in other codes, or to 
examine how galaxy properties depend on merging his- 
tory for example. If such output is requested the mass of 
each node in the tree is recorded, along with the cosmic 
time at which it exists, indices describing relationships 
between parent, child and sibling nodes and, if desired, 
quantities such as the virial radius, virial velocity and 
dark matter scale radius. 

Additionally, Galacticus allows for arbitrary addi- 
tional outputs to be easily implemented. 

2.2.2. Post-processing of Galaxy Properties 

Galacticus is provided with a Perl module that al- 
lows for easy extraction of datasets from a Galacticus 
output file together with a straightforward way to im- 
plement derived properties (i.e. properties computed by 
post-processing the output data). Implementations for 
simple dust-extinction calculations are provided which 
utilize this framework, together with modules that con- 
vert output luminosities into absolute magnitudes in AB 
or Vega systems. Any such derived properties can be 
stored back to the Galacticus output file. 

2.3. Node Evolution Engine 

Galacticus's main task is to evolve galaxies (and 
their associated dark matter halos) through a complex 
merging hierarchy. It begins with a pre-constructed dark 
matter halo merger tree as depicted in the "Stage 1" 
panel of Fig. [T] In this figure, colored circles represent 
nodes in the tree (larger circles implying greater mass 
total mass) and the y-axis represents time increasing up- 
wards (such that t\ > t% > ?3 etc.). Solid lines connect 
main progenitor nodes (typically the most massive) to 



'"Main branch" is defined by starting from the root node of a tree 
and repeatedly stepping back to the most massive progenitor of the 
branch. This does not necessarily select the most massive progenitor 
at a given time. 



their parents, while dotted lines connect other progeni- 
tor nodes (ones that will become subhalos) to their par- 
ents. Each node in the tree is given a unique ID number. 
Galacticus begins with the root node (number 1) and 
checks to see if it can be evolved forward in time. Since 
it has some progenitors, it cannot. Therefore, Galacti- 
cus steps to the primary progenitor (indicated by the 
dashed red arrow) and applies the same condition. This 
eventually leads it to node 4 which has no progenitor 
and so can be evolved forward in time (indicated by the 
dotted green arrow). 

Once node 4 reaches the time at which its parent, 
node 3, is defined, node 4 has effectively become node 

3 and is promoted, replacing node 3. This is a node pro- 
motion event in Galacticus parlance. Since it still has 
no progenitors, Galacticus will continue to evolve node 

4 until it reaches node 2 and is promoted to replace it, as 
indicated in "Stage 2" in Fig.[T] Since node 4 now has 
progenitors it cannot be evolved and Galacticus steps 
back to node 6, which it then evolves forward in time, 
until it is promoted to replace node 5. Following this, 
node 7 is evolved until it reaches node 6. Since node 
7 was not the primary progenitor of node 6, it will be- 
come a subhalo within node 6 (subhalos are indicated 
by dot-dashed circles in Fig. [TJ. This is termed a node 
merger event in Galacticus parlance, and is shown in 
"Stage 3" of Fig. [T] Nodes 6 and 7 are then evolved 
in lockstep until they reach node 4. Galacticus is able 
to handle nested substructure hierarchies, but the cur- 
rent implementation forces there to be a single level of 
subhalos. Therefore, both nodes 6 and 7 now become 
subhalos within node 4. 

Nodes 4, 6 and 7 now evolve in lockstep until "Stage 
4" of Fig. [T] is reached. Here, subhalo node 6 has 
reached the center of its host node 4 (due to the action 
of dynamical friction) and so any galaxies residing in 
nodes 4 and 6 will merge. This is a satellite merger 
event in Galacticus parlance. After this, nodes 4 and 7 
evolve in lockstep until node 4 is promoted to replace 
node 1 as shown in "Stage 5" of Fig. [T] 

Finally, Galacticus descends the remaining branch 
of the merger tree and evolves it forward in time until 
"Stage 6" of Fig.[T]is reached. Since all nodes are now 
at the final time, t\, Galacticus stops evolving them and 
outputs any properties. 

It should be noted that, under this algorithm, at any 
given point in the calculation different branches of the 
tree may have been evolved to different times (as in 
"Stage 5" of Fig. [T] where the left-branch has been 
evolved to t\, but the right-branch has yet to be evolved 
at all). This is acceptable providing that galaxies in dif- 
ferent branches do not interact. If such interactions were 
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Figure I: Description on following page. 



Stage 5 




Figure 1: (cont.) A schematic representation of the evolution of a merger tree within Galacticus. Colored circles represent nodes in the merger 
tree (with larger circles indicating greater mass). Solid circles are isolated halos, while dashed circles represent subhalos inside a larger halo. Time 
progresses from bottom to top (i.e. t\ > t% > t3. , .), Red dashed arrows indicate Galacticus stepping backward through the merger tree to find 
a halo which can be evolved. Green dotted arrows indicate a halo being evolved forward in time. Finally, linked blue circles indicate a merger 
between a subhalo and its host halo and any galaxies that they contain. The evolution of the merger tree is depicted at six key stages, each of which 
is discussed in the text. 
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included it would be necessary to impose a synchronic- 
ity condition on all branches so that no one branch could 
evolve significantly past any other. By design, Galacti- 
cus makes it easy to add in additional time-stepping cri- 
teria such as would be needed to implement such an al- 
gorithm. 

The evolution of nodes forward in time is car- 
ried out by Galacticus's ODE solver, which uses 
GSL's embedded Runge-Kutta-Fehlberg (4, 5) method 
(gsl_odeiv_step_rkf 45). However, any number of 
events are allowed to trigger an "interrupt" which causes 
ODE solving to cease at a specified time. When an in- 
terrupt is generated, a function is specified which will be 
called to handle the interrupt. The function can manipu- 
late the node being evolved as necessary before passing 
it back to the ODE solver. Node promotion, node merg- 
ing and satellite merging events all trigger interrupts. 
Interrupts are also used, for example, to create/destroy 
components as required. 



3. Implemented Components 

In this section we describe the currently implemented 
node components available within Galacticus. We em- 
phasize that a key feature of Galacticus is the abil- 
ity to easily add new components, or replace the cur- 
rent components with alternative implementations. In 
several cases we include "null" implementations which 
effectively switch off a given component. Through- 
out, input parameters are identified using the form 
[inputParameter] . For each component we describe 
the following: 

Properties Any variables (e.g. mass, angular momen- 
tum, etc.) that describe the component; 

Initialization How the component properties are ini- 
tialized prior to any evolution of the merger tree; 

Differential Evolution The equations governing the 
smooth evolution of the component; 

Event Evolution How the component changes in re- 
sponse to one of the following events: 

Node mergers Triggered whenever one dark mat- 
ter halo becomes a subhalo within a large 
halo; 

Satellite merging Triggered whenever a subhalo 
merges with its host halo (or potentially an- 
other subhalo); 



Node promotion Triggered when the main pro- 
genitor of a given node in the merger tree has 
evolved to the time at which its parent node is 
defined, and so must be promoted to become 
that node. 

3.1. ( Supermassive ) Black Hole 

3.1.1. "Null" Implementation 

The null black hole implementation defines the same 
properties as all other black hole implementations, but 
implements dummy functions for all black hole proper- 
ties. It can be used to effectively switch off black holes. 
Of course, this is not safe if any of the other active com- 
ponents expect to get or set black hole properties (or 
if they rely on a sensible implementation of black hole 
evolution). 

3.1.2. "Standard" Implementation 

Properties: The standard black hole implementation 
defines the following properties: 

• The mass of the black hole: M.; 

• The spin of the black hole, j.. 

Initialization: Black holes are not initialized, 
they are created (with a seed mass given by 
[blackHoleSeedMass] and zero spin) as needed. 

Differential Evolution: The mass and spin evolve as: 



M. = (1 - e radiation )M 
= j{M.,j.,M ), 



(1) 
(2) 



where Mq is the rest mass accretion rate, e ra diation is the 
radiative efficiency of the accretion flow feeding the 
black hole and j(M., j.,Mo) is the spin-up function of 
that accretion flow (see §4.3| >. The rest mass accretion 
rate is computed assuming Bondi-Hoyle-Lyttleton 
accretion (see, for example, |Edgar| [2004) > from the 
spheroid gas reservoir (with an assumed temperature of 
[bondiHoyleAccretionTemperatureSpheroid] ) 
enhanced by a factor of 

[bondiHoyleAccretionEnhancementSpheroid] 
and from the host halo (with whatever tem- 
perature the hot halo temperature profile spec- 

S PTTT) 



ifies; 



see 



enhanced by a factor of 
[bondiHoyleAccretionEnhancementHotHalo] . 
The rest mass accretion rate is removed (as a mass sink; 
see j p.4.2| i from the spheroid component. The black 
hole is assumed to cause feedback in two ways: 

Radio-mode Any jet power from the black hole- 
accretion disk system (see §4.3|l is included in the 
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hot halo heating rate providing that the halo is in 
the slow cooling regime (i.e. if the cooling radius 
is smaller than the virial radius; see, for example, 
|Benson and Bowerp OlOa); 



Quasar-mode A mechanical wind luminosity of (Os- 
|trikeretaLl[20T0l ) 



(3) 



where e. wint j — [blackHoleWindEf f iciency] is 
the black hole wind efficiency, is added to the gas 
component of the spheroid (which, presumably, 
will respond with an outflow for example) if and 
only if the wind pressure (at the spheroid charac- 
teristic radius) is less than the typical thermal pres- 
sure in the spheroid gas ( Ciotti etaL]|2009[ ), i.e. 



wind 



1 



,PwindV wind < 



PlSM 

SIcb^ism^ism) 
2ot h 



(4) 



Since Qr 2 p wind V^ nd = L wind where Q is the solid 
angle of the wind flow, this can be rearranged to 
give <pism> > Pwind,criticai where 



Pwind,critical — 



wind 



3Qr 2 y w i nd k B r IS M' 



(5) 



This critical wind density is computed at the char- 
acteristic radius of the spheroid, r sp h e roid, assuming 
Vwind = 10 4 km/s, Tism = 10 4 K and Q. — n, and the 
ISM density is approximated by 



<Pism> 



3A/g asS ph ero id 3 



An 



^"spheroid ' 



(6) 



For numerical ease, the fraction, / w ; nd , of the 
wind luminosity added to the spheroid is adjusted 
smoothly through the p ISM ~ Pwind.criticai region ac- 
cording to 



wind 




if x < o, 
if < x < 1, 
if x > 1, 



(7) 



where x = pisM/Pwind,criticai - 1/2. 

Event Evolution — > Node mergers: None. 

Event Evolution — > Satellite merging: The black 
holes in the two merging galaxies are instantaneously 
merged. Properties are computed using the selected 
black hole binary merger method (see §4.24[ ). 

Event Evolution — > Node promotion: None. 



3.2. Hot Halo 

3.2.1. "Null" Implementation 

The null hot halo implements dummy functions for 
all hot halo properties. It can be used to effectively 
switch off hot halos. Of course, this is not safe if any 
of the other active components expect to get or set hot 
halo properties (or if they rely on a sensible implemen- 
tation of hot halo evolution). 

3.2.2. "Standard" Implementation 

Properties: The standard hot halo implementation 
defines the following properties: 

• The mass of gas which failed to accrete in to the 

hot hal(£] Mfeiied; 

• The mass of gas in the hot halo: M no t; 

• The angular momentum of the gas in the hot halo, 

• The mass(es) of heavy elements in gas in the hot 
halo, M z>h o t ; 

• The mass of gas from outflows in the hot halo: 

^outflowed; 

• The angular momentum of the outflowed gas in the 
hot halo, 7 outflowed ; 

• The mass(es) of heavy elements in outflowed gas, 

-^Z.outflowedi 

and the following pipes: 

Energy Input Energy sent through this pipe is added 
to the hot halo and used to offset the cooling rate 
(see below). 

Cooling Gas The net cooling rate of gas mass (and 
metal content and angular momentum) is sent 
through this pipe. Any component may claim this 
pipe and connect to it, allowing it to receive the 
cooling gas. For example, in the current implemen- 
tation, the exponential disk component (see ^3.3.2 1 
claims and connects to this pipe. 



"By "failed to accrete" we mean any mass which would accrete 
into the halo in the absence of baryonic physics, such as pressure and 
heating. This provides a way to accrete this mass later if, for exam- 
ple, the dark matter halo potential well deepens and becomes more 
effective at accreting mass. 
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Gas Outflow Galactic components that wish to expel 
gas due to an outflow can send that mass (plus 
metals and angular momentum) through this pipe, 
where it will be received into the hot halo compo- 
nent. 

Gas Mass Sink Removes gas (and proportionate 
amounts of angular momentum and elements) 
from the hot gas halo. 

Initialization: At initialization, any nodes with no 
progenitors are assigned a hot halo mass, and failed 
accreted mass as dictated by the baryonic accretion 
method (see § |4.1| > and angular momentum based on the 
accreted mass and the halo spin parameter. 

Differential Evolution: In the standard hot halo 
implementation the hot gas mass and heavy element 
mass(es) evolves as: 



^failed 
M hot 



M; 



Z.hot 



^failed accretion 
^accretion — ^cooling 
^-^outflow,return > 

-m Mz ' hot 

'"cooling . , 

M hot 

^~^Z,outflow,return? 



(8) 
(9) 



(10) 
(11) 



where M accret i on is the rate of growth of the hot com- 
ponent due to accretion from the intergalactic medium 
(IGM) , Mf a ii e d accretion is the rate of failed accretion from 
the IGM (these may include a component due to transfer 
of mass from the failed to hot reservoirs) and M coo ii ng is 
the rate of mass loss from the hot halo due to cooling 



(see §4.5.2 1 minus any heating rate defined as 



M heating = EinpUt/V* 



(12) 



where E is the rate at which energy is being sent through 
the "energy input" pipe and Vyiriai is the virial velocity of 
the haloJ^The angular momentum of the hot gas evolves 
as: 



./h 



J node 



-^node 
outflow.return? 



^cooling ^"cool ^rotate 



(13) 



where M noae and y noQe are defined in 5 3.5 r coo i is the 
cooling radius (see 5 4.5.3| l and V rc .tate is the effective ro- 
tation speed of the halo for its angular momentum (see 



3 The net cooling rate is never allowed to drop below zero. If the 
mass heating rate exceeds the mass cooling rate then the excess energy 
is not used. 



§ |4.6.1 1. For the outflowed components: 



-^outflowed 
-^Z.outnowed 



Moutflow,return ^outflows? (14) 
~^Z,outflow,return +M Zout fl ows , (15) 

(16) 



and: 



' outflowed 



-./, 



outflow,retum J outflows • 



In the above 



,return — Q'outflow return rate 



(17) 



M|M z |7 outflowed 

^"dynamical,halo 



(18) 



where croutflow return rate = [hotHaloOutf lowReturnRate] 
is an input parameter controlling the rate at which gas 
flows from the outflowed to hot reservoirs, and 
M|Mz|7 out fl ows are the net rates of outflow from any 
components in the node. 

Event Evolution — > Node mergers: If the 
[starveSatellites] parameter is true, then any hot 
halo properties of the minor node are added to those of 
the major node and the hot halo component is removed 
from the minor node. Additionally in this case, any 
material outflowed from the the satellite galaxy to its 
hot halo is transferred to the hot halo of the host dark 
matter halo after each timestep. 

Event Evolution — > Satellite merging: If the 
[starveSatellites] parameter is false, then any hot 
halo properties of the satellite node are added to those 
of the host node and the hot halo component is removed 
from the satellite node. 

Event Evolution — > Node promotion: Any hot halo 
properties of the parent node are added to those of the 
node prior to promotion. 

3.3. Galactic Disk 

3.3.1. "Null" Implementation 

The null disk implementation implements dummy 
functions for all disk properties. It can be used to ef- 
fectively switch off disks. Of course, this is not safe if 
any of the other active components expect to get or set 
disk properties (or if they rely on a sensible implemen- 
tation of disk evolution). 

3.3.2. "Exponential" Implementation 

This implementation assumes a disk with an expo- 
nential surface density profile in which stars trace gas. 

Properties: The exponential galactic disk implemen- 
tation defines the following properties: 



The mass of gas in the disk: M^ 



disk.gas > 
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• The masses of elements in the gaseous disk: 

■^Z.disk.gas; 

• The mass of stars in the disk: Mdi s k, s tars; 

• The masses of elements in the stellar disk: 

-^Z,disk,stars> 

• The luminosities (in multiple bands) of the stellar 

disk: -Ldisk.stars; 

• The angular momentum of the disk, /disk; 

• The radial scale length of the disk, 7?disk; 

• The circular velocity of the disk at 7? d isk, Vdisk- 

Initialization: No initialization is performed — disks 
are created as needed. 

Differential Evolution: In the exponential galactic 
disk implementation the gas mass evolves as: 



^^disk.gas — -^cooling -^outflow,disk 
-A^disk.gas/Tbar, 

where the rate of change of stellar mass is 

Mdisk,stars =X P- R- Mdi s k,stars/Tbar, 

with 

* diskjStar formation 

where_Tdisk,star formation is the star formation timescale 



^stars,disk 

(19) 
(20) 

(21) 



(see §417 1, R is the rate of mass recycling from stars 



and Tbar is a bar instability timescale (see §4.7 1. The 



mass removed from the disk by the bar instability mech- 
anism is added to the active spheroid component. El- 
ement abundances (including total metals) evolve ac- 
cording to: 



M; 



Z,disk,gas 



^Z.cooling _ A#Z,outflow,disk 

-M z>starSidisk + y, (22) 



and 



• ^Z,disk,gas • 
™Z,stars,disk - T — K Z 



M, 



(23) 



disk,gas 



where y is the rate of element yield from stars and Rz 
is the rate of element recycling. Recycling rates and 
yields are discussed in §4.12 and §4.18 The angular 
momentum evolves as: 



Isk 



J cooling 

- -^outflow.disk + 
/disk 



M disk , gas + ,V7; 



A^disk.gas + A^disk.stars 



(24) 



The outflow rate, Moutfiw.disk, is computed for the cur- 
rent star formation rate and gas properties by the stellar 
properties subsystem (see § j4.18| l. 

Event Evolution — > Node mergers: None. 

Event Evolution — > Satellite merging: Disks may be 
destroyed (or, potentially, created or otherwise modi- 
fied) as the result of a satellite merging event, as dictated 
by the selected merger remnant mass movement method 
(see j p~97T) . 

Event Evolution — > Node promotion: None. 

3.4. Galactic Spheroid 

3.4.1. "Null" Implementation 

The null spheroid implements dummy functions for 
all spheroid properties. It can be used to effectively 
switch off spheroids. Of course, this is not safe if any of 
the other active components expect to get or set spheroid 
properties (or if they rely on a sensible implementation 
of spheroid evolution). 

3.4.2. "Hernquist" Implementation 

This implementation assumes a Hernquist profile 
( |Hernquist] |1990[ ) for the spheroidal component of a 
galaxy in which stars trace gas. 

Properties: The Hernquist galactic spheroid imple- 
mentation defines the following properties: 

• The mass of gas in the spheroid: M sp h e roid,gas; 

• The masses of elements in the gaseous spheroid: 

-^Z,spheroid,gas > 

• The mass of stars in the spheroid: M sp h e roid,stars; 

• The masses of elements in the stellar spheroid: 

^^Z,spheroid,stars » 

• The luminosities (in multiple bands) of the stellar 
spheroid: L spheroid>stal - s ; 

• The pseudo-angular momentum^] of the spheroid, 

/ spheroid j 

• The radial scale length of the spheroid, r sp heroid; 

• The circular velocity of the spheroid at r sp h er oid, 

^spheroid ■ 

and the following pipes: 



Effectively the angular momentum that the spheroid would have, 
were it rotationally supported rather than pressure supported. 
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Energy Input Energy sent through this pipe is added to 
the gas of the spheroid and will result in an outflow 
(see below). 

Gas Mass Sink Removes gas (and proportionate 
amounts of angular momentum and elements) 
from the spheroid gas. 

Initialization: No initialization is performed — 
spheroids are created as needed. 

Differential Evolution: In the Hernquist galactic 
spheroid implementation the gas mass evolves a^j 



M. 



spheroid, gas 



-M, 



outflow, spheroid 



M 



stars, spheroid; 



where the rate of change of stellar mass is 

^stars, spheroid — ^i* ~ R 

with 

•^spheroid,gas 



(25) 



(26) 



(27) 



spheroid,star formation 



with T SDh eroid.st ar form ation being the star formation 
timescale (see f 4.17 1 and R is the rate of mass recycling 
from stars. Element abundances (including total metals) 
evolve according to: 



^^Z,spheroid,gas — ^^Z.outflow, spheroid ^^Z,stars, spheroid "K) 7 ? (28) 



and 



M: 



Z,stars, spheroid 



m ^^Z.spheroid.gas • 

= i —. Rz 



M, 



(29) 



spheroid,gas 



where y is the rate of element yield from stars and Rz 
is the rate of element recycling. Recycling rates and 
yields are discussed in § |4.12| and §4.18| The angular 
momentum evolves as: 



where Apheroid,energy is set by the 
[spheroidEnergeticOutf lowMassRate] input 
parameter, and £g as ,spheroid is any input energy sent 
through the "Energy Input" pipe. 

Event Evolution — > Node mergers: None. 

Event Evolution — > Satellite merging: Spheroids may 
be created as the result of a satellite merging event, as 
dictated by the selected merger remnant mass move- 
ment method (see §4.9. 1| >. 

Event Evolution — > Node promotion: None. 

3.5. Basic Properties 

Basic properties are the total mass of a node and the 
cosmic time at which it currently exists. 

3.5.1. "Simple" Implementation 

Properties: The simple basic properties implementa- 
tion defines the following properties: 

• The total mass of the node: M no d e ; 
The time at which the node is defined: f no d e ; 

• The time at which the node was last an isolated 
halo (i.e. not a subhalo). 

Initialization: All basic properties are required to be 
initialized by the merger tree construction routine (see 

Differential Evolution: Properties are evolved ac- 
cording to: 



J spheroid — ^^outflow, spheroid 



spheroid 



M s 



spheroid, gas 



• M sph , 



ieroid, stars 



■(30) 



The outflow rate, M utfiow,disk, is computed for the cur- 
rent star formation rate and gas properties by the stel- 
lar properties subsystem (see §4. 18| l, with an additional 
contribution given by 



^^outflow, spheroid — /^spheroid.energy ^2 



gas, spheroid 



(31) 



spheroid 



5 There may be an additional contribution to the mass and angu- 
lar momentum rates of change in the spheroid due to material trans- 
ferred from the disk component via the bar instability mechanism (see 
j |3.3.2( . This is not included here as it is not intrinsic to this specific 
spheroid implementation — it is handled explicitly by the disk compo- 
nent and so applies equally to any spheroid component implementa- 
tion. 



M, 



node 



^node,parenl~^nodc 
^ node, parent — ^node 



if primary progenitoj^ 



otherwise, 



^node — 1 > 



(33) 



where the "parent" subscript indicates a property of the 
parent node in the merger tree. 

Event Evolution — > Node mergers: None. 

Event Evolution — > Satellite merging: None. 

Event Evolution — > Node promotion: M no( j e is up- 
dated to the node mass of the parent prior to promotion. 

3.6. Satellite Node Orbit 

This component tracks the orbital properties of sub- 
halos. 
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3.6.1. "Simple" Implementation 

Properties: The simple satellite orbit implementation 
defines the following properties: 

• The time until the satellite will merge with its host: 

^satellite,merge ■ 

Initialization: None. 

Differential Evolution: Properties are evolved ac- 
cording to: 

^satellite.merge — ~~ 1 ■ (34) 

Event Evolution — > Node mergers: The component is 
created and the time to merging is assigned a value (see 
|4.22.1) . 

Event Evolution — > Satellite merging: None. 
Event Evolution — > Node promotion: Not applicable 
(component only exists for satellite nodes). 

3. 7. Dark Matter Halo Spin 

This component stores and tracks the spin parameters 
of dark matter halos. 

3.7.1. "Null" Implementation 

The null spin component implements dummy func- 
tions for all spin properties. It can be used to effectively 
switch off spins. Of course, this is not safe if any of the 
other active components expect to get or set spin prop- 
erties (or if they rely on a sensible implementation of 
spin evolution). 

3.7.2. "Random" Implementation 

Properties: The random dark matter halo spin imple- 
mentation defines the following properties: 

• The spin parameter of the halo: A. 

Initialization: The spin parameter of each node, if 
not already assigned, is selected at random from a dis- 
tribution of spin parameters. This value is assigned to 
the earliest progenitor of the halo traced along its pri- 
mary branch. The value is then propagated forward 
along the primary branch until the node mass exceeds 
that of the node for which the spin was selected by a 
factor of [randomSpinResetMassFactor] , at which 
point a new spin is selected at random, and the process 
repeated until the end of the branch is reached, in a man- 
ner similar to the algorithm used by |Cole et al.| (|2000). 

Differential Evolution: The spin parameter does not 
evolve. 

Event Evolution — > Node mergers: None 
Event Evolution — > Satellite merging: None. 



Event Evolution — > Node promotion: The spin is up- 
dated to equal that of the parent node. (The two will 
differ only if this is a case where the new halo node was 
sufficiently more massive than the node for which a spin 
was last selected that a new spin value was chosen.) 

3.8. Dark Matter Profile 

This component stores dynamic properties associated 
with dark matter halo density profiles. 

3.8.1. "Null" Implementation 

The null profile implements dummy functions for all 
profile properties. It can be used to effectively switch 
off profiles. Of course, this is not safe if any of the other 
active components expect to get or set profile properties 
(or if they rely on a sensible implementation of profile 
evolution). 

3.8.2. "Scale" Implementation 

Properties: The scale dark matter profile implemen- 
tation defines the following properties: 

• The scale length of the density profile. 

Initialization: The scale length of each 
node, if not already assigned, is assigned us- 
ing the concentration parameter function (see 
j |4.6.2| >, but is not allowed to drop below 
[darkMatterProf ileMinimumConcentration] , 
such that the scale length is equal to the virial radius 
divided by that concentration. The value is propagated 
in both directions along the primary child branch from 
the node. 

Differential Evolution: The scale radius does not 
evolve. 

Event Evolution — > Node mergers: None. 
Event Evolution — > Satellite merging: None. 
Event Evolution — > Node promotion: None. 

4. Specific Physical Implementation 

In this section we describe specific implementations 
of the numerous physical processes that are currently 
available in Galacticus. In each subsection we briefly 
describe what physical process/property is being dis- 
cussed and then describe the specific implementations 
available within Galacticus. 
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4.1. Accretion of Gas into Halos 

The accretion rate of gas from the IGM into a dark 
matter halo is expected to depend on (at least) the rate 
at which that halo's mass is growing, the depth of its 
potential well and the thermodynamical properties of 
the accreting gas. Galacticus implements the follow- 
ing calculations of gas accretion from the IGM. 

Simple Method: Currently the only option, this 
method sets the accretion rate of baryons into a halo to 
be: 







reionization 
Or Z > Zreionization 

otherwise, 



(35) 



where Mhaio is the total rate of growth of the node mass, 
^ionization = [reionizationSuppressionRedshif t] 
is the redshift at which the Universe is reionized and 
Vreionization = [reionizationSuppressionVelocity] 
is the virial velocity below which accretion is sup- 
pressed after reionization. Setting V re ionization to zero 
will effectively switch off the effects of reionization on 
the accretion of baryons. This algorithm attempts to 
offer a simple prescription for the effects of reionization 
and has been explored by multiple authors (e.g. Benson| 
et al. 20021. In particular, |Font et aL ( 20l0| l show 
that it produces results in good agreement with more 
elaborate treatments of reionization. For halos below 
the accretion threshold, any accretion rate that would 
have otherwise occurred is instead placed into the 
"failed" accretion rate. For halos which can accrete, 
and which have some mass in their "failed" reservoir, 
that mass will be added to the regular accretion rate at a 
rate equal to the mass of the "failed" reservoir times the 
specific growth rate of the halo. 

4.2. Background Cosmology 

The background cosmology describes the evolution 
of an isotropic, homogeneous Universe within which 
galaxy formation calculations are carried out. For the 
purposes of Galacticus, the background cosmology is 
used to relate expansion factor/redshift to cosmic time 
and to compute the density of various components (e.g. 
dark matter, dark energy, etc.) at different epochs. 

Matter + Lambda: In this implementation, cosmo- 
logical relations are computed assuming a universe that 
contains only collisionless matter and a cosmological 
constant. 

4.3. Circumnuclear Accretion Disks 

Circumnuclear accretion disks surrounding super- 
massive black holes at the centers of galaxies influence 



the evolution of both the black hole (via accretion rates 
of mass and angular momentum and possibly by extract- 
ing rotational energy from the black hole) and the sur- 
rounding galaxy if they lead to energetic outflows (e.g. 
jets) from the nuclear region. Current implementations 
of accretion disks are as follows and are selected via 
[accretionDisksMethod] . 

Shakura-Sunyaev Geometrically Thin, Radiatively 
Efficient Disks: This implementation assumes that ac- 
cretion disks are always described by a radiatively effi- 
cient, geometrically thin accretion disk as described by 
Shakura and Sunyaev]( |1973| >. The radiative efficiency 
of the flow is computed assuming that material falls into 
the black hole without further energy loss from the in- 
nermost stable circular orbit (IS CO) , while the spin- 
up rate of the black hole is computed assuming that the 
material enters the black hole with the specific angular 
momentum of the ISCO (i.e. there are no torques on the 
material once it begins to fall in from the ISCO;[Bardeen| 



1970]>. For these thin disks, jet power is computed, us- 



ing the expressions from Meier ( 2001} his equations 4 
and 5). 

Advection Dominated, Geometrically Thick, Radia- 
tively Inefficient Flows (ADAFs): This implementa- 
tion assumes that accretion is via an advection dom- 
inated accretion flow (Naray an and Yi| |1994| > which 
is radiatively inefficient and geometrically thick. The 
radiative efficiency of the flow, which will be zero 
for a pure ADAF, can be set via the input parame- 
ter [adaf RadiativeEf f iciency] . The spin up rate 
of the black hole and the jet power produced as ma- 
terial accretes into the black hole are computed using 
the method of Benson a nd Babul] ( |2009[ ). The energy 
of the accreted material can be set equal to the energy 
at infinity (as expected for a pure ADAF) or the en- 
ergy at the ISCO by use of the [adaf EnergyOption] 
parameter. The ADAF structure is controlled by 
the adiabatic index, y, and viscosity parameter, a, 
which are specified via the [adaf Adiabaticlndex] 
and [adaf ViscosityOption] input parameters re- 
spectively, [adaf ViscosityOption] may be set to 
"f it", in which case the fitting function for a as a func- 
tion of black hole spin from Benson a nd Babul] ( |2009| l 
will be used. 

"Switched" Disks: This method allows for accretion 
disks to switch between radiatively efficient (Shakura- 
Sunyaev) and inefficient (ADAF) modes. Which mode 
is used is determine by the accretion rate onto the disk: 

• Radiatively efficient accretion if 
A^/^Eddington > [accretionRateThinDiskMinimum] 
and 
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A^/^Eddington < [accretionRateThinDiskMaximum] ; 
• Radiatively inefficient accretion otherwise. 

4.4. Cold Dark Matter Structure Formation 

A variety of functions are used to describe struc- 
ture formation in cold dark matter dominated universes. 
These are described below. 

4.4.1. Primordial Power Spectrum 

The functional form of the primordial dark matter 
power spectrum. The power spectrum is computed from 
the specified primordial power spectrum and the trans- 
fer function (see §4.4.2| > and normalized to a value of cr 8 
specified by [sigma_8]. 

(Running) Power Law Spectrum: This method imple- 
ments a primordial power spectrum of the form: 



P(k) oc k" M , 



where 



n eS (k) = n s + 



dn k 
dink k m f ' 



(36) 



(37) 



where n s — [power Spectrumlndex] is 

the power spectrum index at wavenumber 
k K f — [powerSpectrumRef erenceWavenumber] 
and dn /din k = [powerSpectrumRunning] describes 
the running of this index with wavenumber. 

4.4.2. Transfer Function 

The functional form of the cold dark matter transfer 
function is selected by [transf erFunctionMethod] . 
The power spectrum is computed from the specified 
transfer function and the primordial power spectrum 
(see §4.4. 1[ ) and normalized to a value of cr% specified 
by [sigma_8] . 

BBKS: This method uses the fitting function of 
Bardeen et al. ( 1986| l to compute the CDM transfer 
function. 

Eisenstein & Hu: This method uses the fitting 
function of Eisens teTn and H"u| ( |1999[ l to compute 
the CDM transfer function. It requires that the ef- 
fective number of neutrino species be specified via 
the [ef f ectiveNumberNeutrinos] parameter and 
summed mass of all neutrino species (in eV) be spec- 
ified via the [summedNeutrinoMasses] parameter. 

CMBFast: This method uses the CMBFast code to 
compute the CDM transfer function. It requires that the 
mass fraction of helium in the early Universe be speci- 
fied via the [Y_He] parameter. CMBFast will be down- 
loaded and run if the transfer function needs to be com- 
puted. The transfer function will then be stored in a file 




Wavenumber [Mpc" ] 



Figure 2: The transfer function, T(k), computed automatically by 
Galacticus by downloading, compiling and running CMBFast. The 
transfer function is computed for a cosmological model with £Im = 
0.25, n A = 0.75, H b = 0.045 and H = 75 km/s/Mpc. 



for future reference. An example of a transfer function 
computed in this way is shown in Fig. [2] 

File: This method reads a tabulated transfer func- 
tion from an XML file, interpolating between tabulated 
points. 

4.4.3. Linear Growth Function 

The function describing the amplitude of linear per- 
turbations. 

Simple: This method calculates the growth of lin- 
ear perturbations using standard perturbation theory in 
a Universe consisting of collisionless matter and a cos- 
mological constant. 

4.4.4. Critical Overdensity 

The method used to compute the critical linear over- 
density at which overdense regions virialize. 

Spherical Collapse (Matter + Cosmological Con- 
stant): This method calculates critical overdensity using 
a spherical top-hat collapse model assuming a Universe 
which contains collisionless matter and a cosmological 
constant (see, for example, Percival|2005 1. 



4.4.5. Virial Density Contrast 

The method used to compute the mean density con- 
trast of virialized dark matter halos is specified via 
[virialDensityContrastMethod] . 

Bryan & Norman Fitting Function: This method cal- 
culates virial density contrast using the fitting functions 
given by Bryan and Norman ( 1998| l. As such, it is valid 
only for Sl\ = or Q.m + = 1 cosmologies and will 
abort on other cosmologies. 
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Spherical Collapse (Matter + Cosmological Con- 
stant): This method calculates virial density contrast us- 
ing a spherical top-hat collapse model assuming a Uni- 
verse which contains collisionless matter and a cosmo- 
logical constant (see, for example, Percival|2005| l. 



4.4.6. Halo Mass Function 

The dark matter halo mass function (i.e. the num- 
ber of halos per unit volume per unit mass inter- 
val). Which mass function to use is specified by 
[haloMassFunctionMethod] . 

Press-Schechter: This method uses the functional 
form proposed by |Press and~S chechter ( |1974[ ) to com- 
pute the halo mass function. 

Sheth-Tormen: This method uses the functional form 



proposed by Sheth et al. (2001 1 to compute the halo 
mass function. 

Tinker: This method uses the functional form pro- 
posed by Tink er et al.| ( |2010) to compute the halo mass 
function. The mass function is computed at the appro- 



priate virial overdensity (see 54.4.5 I. 



4.5. Cooling of Gas Inside Halos 

The cooling of gas within dark matter halos is con- 
trolled by a number of different algorithms which will 
be described below. 

4.5.1. Cooling Function 

The cooling function of gas, A(p, T, Z), where 
p is gas density, T is temperature and Z is 
a vector of elemental abundances, is selected by 
[coolingFunctionMethod] . Multiple such methods 
may be specified and are cumulative (i.e. the net cooling 
function is the sum over all specified cooling functions). 

Atomic Collisional Ionization Equilibrium Using 
Cloudy: This method computes the cooling function us- 
ing the Cloudy code and under the assumption of colli- 
sional ionization equilibrium with no molecular contri- 
bution. Abundances are Solar, except for zero metallic- 
ity calculations which use Cloudy's "primordial" metal- 
licity. The helium abundance for non-zero metallicity is 
scaled linearly with metallicity between primordial and 
Solar values. The Cloudy code will be downloaded and 
run to compute the cooling function as needed, which 
will then be stored for future use. As this process is 
slow, a precomputed table is provided with Galacticus. 
If metallicities outside the range tabulated in this file 
are required it will be regenerated with an appropriate 
range. Fig. [3] shows an example of a cooling function 
computed automatically by Galacticus using Cloudy, 

Collisional Ionization Equilibrium From File: This 
method the cooling function is read from a file. The 




n 



W 



Temperature [K] 



Figure 3: The cooling function for atomic gas in collisional ioniza- 
tion equilibrium computed using Cloudy 08.00 (automatically down- 
loaded, compiled and run by Galacticus) as a function of gas temper- 
ature. The color scale indicates the gas metallicity. Solar abundance 
ratios are assumed, except for zero metallicity calculations which use 
Cloudy's "primordial" metallicity. 



cooling function is assumed to be computed under con- 
ditions of collisional ionization equilibrium and there- 
fore to scale as p 2 . 

CMB Compton Cooling: This method computes the 
cooling function due to Compton scattering off of cos- 
mic microwave background (CMB) photons: 



4cr x akBn e 4 
A = ^cmb (-* _ Tcmb) , 



(38) 



where cr T is the Thompson cross-section, a is the ra- 
diation constant, ke is Boltzmann's constant, n e is the 
number density of electrons, m e is the electron mass, 
c is the speed of light, 7cmb is the CMB temperature at 
the current cosmic epoch and T is the temperature of the 
gas. The electron density is computed from the selected 
ionization state method (see § |4.13| l. 

4.5.2. Cooling Rate 

The algorithm used to compute the rate at which gas 
drops out of the hot halo due to cooling. 

White & Frenk: This method computes the cooling 
rate using the expression given by |White and Frenk| 
( [T99T) , namely 



M coo i = 47ir cool p(r coo i)r 

cool 



(39) 



where r coo i is the cooling radius (see 54.5.3 1 in the hot 
halo and p(r) is the density profile of the hot halo (see 
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4.5.3. Cooling Radius 

The algorithm used to compute the cooling radius in 
hot halo gas. 

Simple: This method computes the cooling radius by 
seeking the radius at which the time available for cool- 



ing (see 54.5.51 equals the cooling time (see 5 4.5.4 1. 
The growth rate is determined consistently based on the 
slope of the density profile, the density dependence of 
the cooling function and the rate at which the time avail- 
able for cooling is increasing. This method assumes that 
the cooling time is a monotonic function of radius. 



4.5.4. Cooling Time 

The algorithm used to compute the time taken for gas 
to cool (i.e. the cooling time). 

Simple: This method assumes that the cooling time is 
simply 



fcool 



N k B Tn u 
2 A 



(40) 



where N = [coolingTimeSimpleDegreesOf Freedom] 
is the number of degrees of freedom in the cooling 
gas which has temperature T and total particle number 
density n tot and A is the cooling function (see §4.5. 

4.5.5. Time Available for Cooling 

The method used to determine the time available for 
cooling (i.e. the time for which gas in a halo has been 
able to cool). 

White & Frenk: This method assumes that the time 
available for cooling is equal to 

^available = exp [/In funiveise + (1 - /) In fdynamical] > (41) 

where / = [coolingTimeAvailableAgeFactor] is 
an interpolating factor, fjjniverse is the age of the Universe 
and fdynamical is the dynamical time in the halo. The orig- 



inal White and Frenk ( 1991[ ) algorithm corresponds to 
/=!• 

4.6. Dark Matter Halos 

Several algorithms are used to implement dark matter 
halos. These are described below. 

4.6.1. Density Profile 

The method used to compute density pro- 
files of dark matter halos is selected via the 
[darkMatterProf ileMethod] parameter. 

Isothermal: Under this method the density profile is 
given by: 



PnodeM °C r 



(42) 



normalized such that the total mass of the node is en- 
closed with the virial radius. 

NFW: Under this method the Navarro-Frenk-White 



(NFW) density profile (Navarro et al. 1997 lis used 



PnodeM K 



1 



(43) 



normalized such that the total mass of the node is en- 
closed with the virial radius and with scale length r s = 
'"vMai/c where c is the halo concentration (see §4.6.2| >. 

4.6.2. Dark Matter Density Profile Concentration 
The method used to compute the concentration of 

dark matter halos. 

Gao2008: Under this method the concentration is 

computed using a fitting function from |Gao et IE] 

( |20Dg) : 



log 10 c = A log 10 M ha io + B. 



(44) 



The parameters are a function of expansion factor, a. 
We use the following fits to the |Gao et al.| ( |2008] l results: 



A = -0.14 exp 
B = 2.646 exp 



log 10 a + 0.05 



log lo a 
0.50 



0.35 

2' 



(45) 
(46) 



4.6.3. Spin Parameter Distribution 

The method used to compute the distribution of 
dark matter halo spin parameters is selected via the 
[haloSpinDistributionMethod] parameter. 

Lognormal: Under this method the spin is drawn 
from a lognormal distribution. 

Bett2007: Under this method the spin is drawn from 
the distribution found by |Bett et al.| ( |2007| >. The A 
and a parameter of Bett et al.'s distribution are set 
by the [spinDistributionBett2007LambdaO] and 
[spinDistributionBett2007Alpha] input parame- 
ters. 

4. 7. Disk Stability/Bar Formation 

The method uses to compute the bar insta- 
bility timescale for galactic disks is selected by 
[barlnstabilityMethod] . 

Efstathiou, Lake & Negroponte: This method uses the 
stability criterion of Efstathio u et al.| ( [T9 82) to estimate 
when disks are unstable to bar formation: 



y peak 



VGMdiskMisk 



< 6c, 



(47) 
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for stability, where y pea k is the peak velocity in the ro- 
tation curve (computed here assuming an isolated ex- 
ponential disk), Mdi s k is the mass of the disk and 
is its scale length (assuming an exponential disk). The 
value of e c is linearly interpolated in the disk gas frac- 
tion between values for purely stellar and gaseous disks 
as specified by [stabilityThresholdStellar] and 
[stabilityThresholdGaseous] respectively. For 
disks which are judged to be unstable, the timescale for 
bar formation is taken to be 



^bar — ^disk 



(48) 



where ei so is the value of e for an isolated disk and fdisk is 
the disk dynamical time, defined as r/V (where V is the 
circular velocity) at one scale length. This form gives 
an infinite timescale at the stability threshold, reducing 
to a dynamical time for highly unstable disks. 

4.8. Galactic Structure 

The algorithm to be used when solving 
for galactic structure (specifically, finding 
radii of galactic components) is specified via 
[galacticStructureRadiusSolverMethod] . 

Simple: This method determines the sizes of galactic 
components by assuming that their self-gravity is neg- 
ligible (i.e. that the gravitational potential well is dom- 
inated by dark matter) and that, therefore, baryons do 
not modify the dark matter density profile. The radius 
of a given component is then found by solving 



j= VGM DM (r)r, 



(49) 



where j is the specific angular momentum of the com- 
ponent (at whatever point in the profile is to be solved 
for), r is radius and M(r) is the mass of dark matter 
within radius r. 

Adiabatic: This method takes into account the 
baryonic self-gravity of all galactic components when 
solving for structure and additionally accounts for 
back reaction of the baryons on the dark matter 
density profile using the adiabatic contraction al- 
gorithm of |Gnedin et al.| ( |2004| l. The parame- 
ters A and a> of that model are specified via in- 
put parameters [adiabaticContractionGnedinA] 
and [adiabaticContractionGnedinOmega] respec- 
tively. Solution proceeds via an iterative procedure to 
find equilibrium radii for all galactic components in a 
consistently contracted halo. The method used follows 
that described by |Benson and B ower (2 010b] >. 



4.9. Galaxy Merging 

The process of merging two galaxies currently in- 
volves two algorithms: one which decides how the 
merger causes mass components from both galaxies to 
move and one which determines the size of the remnant 
galaxy. 

4.9.1. Mass Movements 

The movement of mass elements in the merging 
galaxies. In the following, M\ and M2 are the baryonic 
masses of the satellite and central galaxies respectively 
that are about to merge. 

Simple: This method implements mass movements 
according to: 

• If Mi > f m - d j or M2 then all mass from both satellite 
and central galaxies moves to the spheroid compo- 
nent of the central galaxy; 

• Otherwise: Gas from the satellite moves to 
the component of the central specified by the 
[minorMergerGasMovesTo] parameter (either 
"disk" or "spheroid"), stars from the satellite 
move to the spheroid of the central and mass in the 
central does not move. 

Here, / ma j or = [majorMergerMassRatio] is the mass 
ratio above which a merger is considered to be "major". 

4.9.2. Remnant Sizes 

The method used to calculate the sizes of merger rem- 
nants. 

Null: This is a null method which does nothing at all. 
It is useful, for example, when running Galacticus to 
study dark matter only (i.e. when no galaxy properties 
are computed). 

Cole et al. (2000): This method uses the algorithm of 



Cole et al. ( 2000 1 to compute merger remnant spheroid 



sizes. Specifically 
(Mi + M 2 ) 2 



Mf 



M 



2 /orbit M V M 2 



(50) 



rnew n r 2 c n + r 2 

where Mj and M2 are the baryonic masses of the merg- 
ing galaxies and r\ and r2 are their half mass radii, r new 
is the half mass radius of the spheroidal component of 
the remnant galaxy and c is a constant which depends on 
the distribution of the mass. For a Hernquist spheroid 
c = 0.40 can be found by numerical integration while 
for a exponential disk c = 0.49. For simplicity a value 
of c = 0.5 is adopted for all components. The param- 
eter / or bit = [mergerRemnantSizeOrbitalEnergy] 
depends on the orbital parameters of the galaxy pair. For 
example, a value of /orbit = 1 corresponds to point mass 
galaxies in circular orbits about their center of mass. 
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4.10. Hot Halo Density Profile 

The hot halo density profile is used, for example, in 
calculations of cooling rates. 

Cored Isothermal: This method adopts a spherically 
symmetric cored-isothermal density profile for the hot 
halo. Specifically, 



Phot haloM <* [r 2 + ^core] 



(51) 



where the core radius, r oore , is set to be a 
fixed fraction of the virial radius, that frac- 
tion being given by the input parameter 
[isothermalCoreRadiusOverVirialRadius] . 
The profile is normalized such that the current mass in 
the hot gas profile is contained within the virial radius. 

4.11. Hot Halo Temperature Profile 

The hot halo temperature profile is used, for example, 
in calculations of cooling rates. 

Virial Temperature: This method assumes an isother- 
mal halo with a temperature equal to the virial tempera- 
ture of the halo. 

4.72. Initial Mass Functions 

The stellar initial mass function (IMF) subsystem 
within Galacticus supports multiple IMFs and extensi- 
ble algorithms to select which IMF to use based on the 
physical conditions of star formation. 

4.12.1. Initial Mass Function Selection 

The method to use for selecting which IMF to use. 

Fixed: This method uses a fixed IMF irrespective of 
physical conditions. The IMF to use is specified by the 
[imf SelectionFixed] parameter (e.g. setting this pa- 
rameter to Salpeter selects the Salpeter IMF). 

4.12.2. Initial Mass Functions 

A variety of different IMFs are available. Each 
IMF supplies a recycled fraction and metal yield 
for use in the instantaneous recycling approx- 
imation. These can be set via the parameters 
imf {imf Name}RecycledInstantaneous and 
imf {imf Name}YieldInstantaneous where 
{imf Name} is the name of the IMF. 

Chabrier: The Chabrier IMF is defined by 
HChabrier||200l) : 



<p(M) oc < 



AT'e 



[log 10 (f)/ C r c ] 2 /2 



M 



-2.3 



for0.1M o 

< M < \M B 
for 1M 

< M < 125M 
otherwise, 



where cr c = 0.69 and M c = 0.08M Q . 

Kennicutt: The Kennicutt IMF is defined by (Ken- 
[mcuttl[T983l l: 
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(53) 



Kroupa: The Kroupa IMF is defined by (Kroupa 
12001) : 
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(54) 



Miller-Scalo: The Miller-Scalo IMF is defined by 
dMillerand Scalol [19791 : 



4>{M) ■ 



' M 


-1.25 


forO.lOM 


< M < l.OOM 


M 


-2.00 


for 1.00M o 


< M < 2.00M G 


M~ 


-2.30 


for 2.00M Q 


< M < 1O.OM (55) 


M~ 


-3.30 


for 10.0M G 


< M < 125M 


, 




otherwise. 





Salpeter: The Salpeter IMF is defined by ( |Salpeter| 
[1955): 
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(56) 



Scalo: The Scalo IMF is defined by (Scalo 19861: 



4>{M) ■ 



( M +im 


forO.lOM 


< M < O.18M 


M~ 


-1.01 


forO.18M 


< M < O.42M 


M' 


-2.75 


for O.42M 


< M < O.62M 


M' 


-2.08 


for O.62M 


< M < 1.18M (57) 


M~ 


-3.50 


for 1.1 8M 


< M < 3.50M o 


M~ 


-2.63 


for 3.5OM 


< M < 125M 


, 




otherwise. 





4.13. Ionization State 

The ionization state of gas, including, for example, 
the electron density, as a function of the gas density, 
composition and temperature. 

Atomic Collisional Ionization Equilibrium Using 
Cloudy: This method computes the ionization state us- 
ing the Cloudy code and under the assumption of colli- 
sional ionization equilibrium with no molecular contri- 
bution. Abundances are Solar, except for zero metallic- 
(52) ity calculations which use Cloudy's "primordial" metal- 
licity. The helium abundance for non-zero metallicity is 
scaled linearly with metallicity between primordial and 
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Solar values. The Cloudy code will be downloaded and 
run to compute the cooling function as needed, which 
will then be stored for future use. As this process is 
slow, a precomputed table is provided with Galacticus. 
If metallicities outside the range tabulated in this file 
are required it will be regenerated with an appropriate 
range. 

Collisional Ionization Equilibrium From File: In this 
method the ionization state is read from a file. The ion- 
ization state is assumed to be computed under condi- 
tions of collisional ionization equilibrium and therefore 
densities of all species scale as the total density, p. 

4.14. Merger Tree Construction 

Merger trees are "constructed' using a method spec- 
ified via [mergerTreeConstructMethod] . 

Read From File: This method reads merger tree struc- 
tures from an HDF5 file. 

Build: This method first creates a distribution of tree 
root halo masses at a specified final redshift and then 
builds a merger tree using the selected build algorithm 
(see §4.16[ ). The root halo masses are selected to lie 
within a user specified range, with a specified average 
number of trees per decade of halo mass. The distribu- 
tion of halo masses is such that the mass of the 2 th halo 
is 

M hal0ii = exp [ln(M hal0;min ) 

+ In (M ha i 0>max /M ha i 0>min ) x] + °] . (58) 

Here, jc, is a number between and 1 and a is an in- 
put parameter that controls the relative number of low 
and high mass tree produced. The distribution of x is 
controlled by an input parameter with options: 

uniform x is distributed uniformly between and 1; 

quasi x is distributed using a quasi-random sequence. 

4.15. Merger Tree Branching 

The method to be used for computing branching 
probabilities in merger trees when trees are constructed 
using Monte Carlo techniques. 

Modified Press-Scheduler: This method uses 
the algorithm of Park inson et al.| ( |2008| ) to com- 
pute branching ratios. The parameters Go, j\ 
and ji of their algorithm are specified by the in- 
put parameters [modif iedPressSchechterGO] , 
[modif iedPressSchechterGammal] and 



[modif iedPressSchechterGamma2] re- 
spectively. Additionally, the parameter 
[modif iedPressSchechterFirstOrderAccuracy] 
limits the step in the critical linear theory over- 
density for collapse, <5 cr j t , so that it never exceeds 
[modif iedPressSchechterFirstOrderAccuracy] 
times ^2[o- 2 (M 2 /2) - cr 2 (M 2 )], where M 2 is the mass 
of the halo being considered for branching and cr(M) 
is the CDM mass variance computed by filtering the 
power spectrum using top-hat spheres. This ensures 
that the first order expansion of the merging rate that 
is assumed is accurate. Progenitor mass functions 
computed by this branching algorithm when used in the 
|Cole et aT] ( |2000| l merger tree building algorithm are 
shown in Fig. [4] where they are compared with equiv- 
alent mass functions measured from the Millennium 
Simulation (Springel et al. 2005 1. Clearly there is 



By "construct" we mean any process of creating a representation 
of a merger tree within Galacticus. 



excellent agreement with the N-body results. 

4.76. Merger Tree Building 

The method to be used for building merger trees. 

Cole et al. (2000) Algorithm: This method 
uses the algorithm described by Col e et aT] ([2000), 
with a branching probability method selected via the 
treeBranchingMethod parameter. This action of this 
algorithm is controlled by the following parameters: 

[mergerTreeBuildCole2000MergeProbability] 
The maximum probability for a binary merger 
allowed in a single timestep. This allows the prob- 
ability to be kept small, such the the probability 
for multiple mergers within a single timestep is 
small. 

[mergerTreeBuildCole2000AccretionLimit] 

The maximum fractional change in mass due to 
sub-resolution accretion allowed in any given 
timestep when building the tree. 

[mergerTreeBuildCole2000MassResolution] 

The minimum halo mass that the algorithm will 
follow. Mass accretion below this scale is treated 
as smooth accretion and branches are truncated 
once they fall below this mass. 

Smooth Accretion: This method builds a branchless 
merger tree with a smooth accretion history using the 
fitting formula of Wechsler et al. (2002). The tree is de- 
fined by a final mass at a specified redshift and is con- 
tinued back in time by decreasing the halo mass by a 
specified factor at each new node until a specified mass 
resolution is reached. The fitting formula of |Wechsler| 
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□ 

Figure 4: Progenitor mass functions at redshifts z = 0.5, 1, 2 and 4 (bottom to top) for halos of mass 10 12±0151 , 10 13 - 5±0 151 and 10 I5±alsI A _1 Af o 
(left to right) at z = are shown. Here, M\ is the mass of the progenitor halo and M2 the mass of the z = halo. Green lines were measured from 
the Millennium Simulation ( Springel et al. 2005 1 by Parkinson et al. (2008 1, while red lines are computed using Galacticus's merger tree building 
routines (with the[Parkinson et al. (2008 1 branching algorithm and the Cole et al. (2000) tree building algorithm). 
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et al. ( 2002| ) has one free parameter, the formation red- 
shift. The formation redshift can either be computed au- 
tomatically using the method of Bullock et al. (2001), or 
specified by an input parameter. 

4.17. Star Formation Timescales 

The methods for computing star formation timescales 
in disks and spheroids. 

Power Law: This method (as used by |Cole et aL| 
(2000) for example) computes the star formation 



populations are computed, for any given IMF, from stel- 
lar evolution models. The rates of change are then: 



M+ = 



-f 
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cf>(t')R(t-t';Z fuel [t'])dt', (61) 
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timescale to be: 



' ★ — ^"dynamical 



V 



200km/s 



(63) 



(59) 



where = [starFormation-[Component}Ef f iciency] 
and a* — [starFormation{Component}VelocityExponent] 
are input parameters (with {Component} being a place- 
holder for Disk or Spheroid), Td ynam i ca i = r/V 
is the dynamical timescale of the component and 
r and V are the characteristic radius and velocity 
respectively of the component. The timescale is not 
allowed to fall below a minimum value specified by 
[starFormation{Component}MinimumTimescale] . 



4.18. Stellar Population Properties 

The algorithm for determining stellar population 
properties — essentially the rates of change of stellar and 
gas mass and abundances given a star formation rate 
and fuel abundances (and perhaps a historical record 
of star formation in the component)-is specified by 
[stellarPopulationPropertiesMethod] . 

Instantaneous: This method uses the instantaneous 
recycling approximation. Specifically, given a star for- 
mation rate <p, this method assumes a rate of increase of 
stellar mass of M* = (1 - R)(f>, and a corresponding rate 
of decrease in fuel mass. The rate of change of the metal 
content of stars follows from the fuel metallicity, while 
that of the fuel changes according to 



M 



fud,Z 



-(1 -R)Z {ud <p + y<p. 



(60) 



In the above R is the instantaneous recycled frac- 
tion and y is the yield, both of which are sup- 
plied by the IMF subsystem (see 54.12). The rate 



of energy input from the stellar population is com- 
puted assuming that the canonical amount of en- 
ergy from a single stellar population (as defined by 
the f eedbackEnergylnputAtlnf inityCanonical) 
is input instantaneously. 

Non-instantaneous: This method assumes fully non- 
instantaneous recycling and metal enrichment. Recy- 
cling and metal production rates from simple stellar 
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where R(t;Z) and p(t;Z) are the recycling and metal 
yield rates respectively from a stellar population of age 
t and metallicity Z. The energy input rate is computed 
self-consistently from the star formation history. 



4.79. Stellar Population Spectra 

Stellar population spectra are used to construct inte- 
grated spectra of galaxies. 

Conroy, White & Gunn: This method uses v2.1 of 
the FSPS code of Conroy et al. (2009) to compute stel- 
lar spectra. If necessary, the FSPS code will be down- 
loaded, patchecj^and compiled and run to generate spec- 
tra. These tabulations are then stored to file for later 
retrieval. 

File: This method reads stellar population spectra 
from an HDF5 file. 

4.20. Stellar Population Spectra Post-processing 

Stellar population spectra are post-processed (to han- 
dle, for example, absorption by the IGM) by any number 
of the following algorithms. 

Meiksin (2006) IGM Attenuation: This method post- 
processes spectra through absorption by the IGM using 
the results of |Meiksin| ( |2006| ). 

Madau (1995) IGM Attenuation: This method post- 
processes spectra through absorption by the IGM using 
the results of |Madau| ( |1995| ). 

Null Method: This method performs no post- 
processing. 



'Patching is not to fix any bugs in FSPS, but merely to insert code 
for reading a tabulated IMF output by Galacticus. 
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4.21. Stellar Astrophysics 

Various properties related to stellar astrophysics are 
required by Galacticus. The following documents their 
implementation. 

4.21.1. Basics 

This subset of properties include recycled mass, 
metal yield and lifetime. 

File: This method reads properties of individual stars 
of different initial mass and metallicity from an XML 
file and interpolates in them. The stars can be irregularly 
spaced in the plane of initial mass and metallicity. 



any Population Ill-specific supernovae energy plus the 
integrated energy input from stellar winds. The number 
of Type II SNe is computed automatically from the IMF 
and a specified minimum mass required for a Type II 
supernova. 

4.22. Substructure and Merging 

Substructures and merging of nodes/substructures is 
controlled by several algorithms which are described 
below: 



4.27.2. Stellar Winds 

Energy input to the ISM from stellar winds is used in 
calculations of feedback efficiency. 

Leitherer et al. (1992): This method uses the fit- 
ting formulae of Leither er et al.| ( |1992| l to compute stel- 
lar wind energy input from the luminosity and effective 
temperature of a star (see §4.21.3[ ). 

4.21.3. Stellar Tracks 

The method used to compute stellar evolutionary 
tracks. 

File: This method luminosities and effective temper- 
atures of stars are computed from a tabulated set of stel- 
lar tracks. Galacticus is supplied with a suitable file 
constructed from the tracks of Ber teili et al.| ( |2008] l and 
|Bertelli et al-1p009) . 

4.21.4. Supernovae Type la 

Properties of Type la supernovae, including the cu- 
mulative number occurring and metal yield. 

Nagashima et al. (2005) Prescription: This method 
uses the prescriptions from Nagashima et al.| ( |2005[ ) to 
compute the numbers and yields of Type la supernovae. 

4.21.5. Population III Supernovae 

Properties of Population III specific supernovae, in 
particular the energy released. 

Heger & Woosley (2002): This method computes the 
energies of pair instability supernovae from the results 



of Heger and Woosley (2002). 



4.21.6. Stellar Feedback 

Aspects of stellar feedback, such as the energy input 
to the surrounding ISM. 

Standard: This method assumes that the cumulative 
energy input from a stellar population is equal to the 
total number of (Type II and Type la) supernovae mul- 
tiplied by a specified energy per supernovae (SNe) plus 



4.22.1. Merging Timescales 

The method used to compute merging 
timescales of substructures is specified via the 
[satelliteMergingMethod] parameter. 

Dynamical Friction: Lacey & Cole: This method 
computes merging timescales using the dynami- 
cal friction calculation of [Lacey and Cole| ( |1993| l. 
Timescales are multiplied by the value of the 
[mergingTimescaleMultiplier] input parameter. 

Dynamical Friction: Jiang (2008): This method 
computes merging timescales using the dynamical fric- 
tion calibration of |Jiang et al.| (2008 ). 

Dynamical Friction: Boylan-Kolchin (2008): This 
method computes merging timescales using the dynam- 
ical friction calibration of Boylan-Kolchin et al. ( 2008 ). 



4.22.2. Virial Orbits 

The algorithm to be used to determine orbital param- 
eters of substructures when they first enter the virial ra- 
dius of their host. 

Benson (2005): This method selects orbital param- 



eters randomly from the distribution given by [Benson 
( |2005) . 



4.22.3. Node Merging 

The algorithm to be used to process nodes when they 
become substructures. 

Single Level Hierarchy: This method maintains a sin- 
gle level hierarchy of substructure, i.e. it tracks only 
substructures, not sub-substructures or deeper levels. 
When a node first becomes a satellite it is appended to 
the list of satellites associated with its host halo. If the 
node contains its own satellites they will be detached 
from the node and appended to the list of satellites of 
the new host (and assigned new merging times; see 
§4.22.1> . 
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4.23. Supernovae Feedback Models 

The supernovae feedback driven outflow rate for 
disks and spheroids. 

Power Law: This method assumes an outflow rate of: 



M, 



outflo 



Voutflo 



V V 



p • < 65 ) 

-'-'canonical 

where Voutflow — [{component}Outf lowVelocity] 
and ffoutflow = [{component}Outf lowExponent] are 
input parameters, V is the characteristic velocity of the 
component, E is the rate of energy input from stellar 
populations and ^canonical is the total energy input by a 
canonical stellar population normalized to 1M after in- 
finite time. 

4.24. Supennassive Black Holes Binary Mergers 

The method to be used for computing the effects of 
binary mergers of supermassive black holes. 

Rezzolla et al. (2008): This method uses the fitting 
function of Rezzolla et al. (2008} to compute the spin 
of the black hole resulting from a binary merger. The 
mass of the resulting black hole is assumed to equal the 
sum of the mass of the initial black holes (i.e. there is 
negligible energy loss through gravitational waves). 

5. Example Calculations 

In this section we show results from an example 
Galacticus model. The parameters of this model are 
listed in Table Q] It should be noted that the values 
of many of these parameters were chosen on the ba- 
sis of previous experience with similar models — we 
have not performed an extensive search of the available 
parameter space, either manually or in an automated 
manner. As such, we do not claim that this example 
model represents the best match to observational data 
that can be obtained with the standard implementation 
of Galacticus — indeed it is certainly not so, although it 
does achieve reasonably good matches to many z ~ 
observational datasets. Instead, this example model 
merely serves to illustrate the type of predictions which 
can be extracted from Galacticus. A detailed search of 
the available parameter space will be presented in a fu- 
ture paper. 

5.1. Numerical Results 

The numerical solution of ODEs in Galacticus is 
controlled by various parameters which affect the ac- 
curacy of solution. In this subsection we explore the 
numerical convergence of Galacticus with respect to 
these results. Additionally, we examine the time taken 
to evolve merger trees in the example model. 



5.7.7. Convergence 

To test the numerical convergence in Galacticus we 
run a single 1O 12 M merger tree to z — as our base 
model. We repeat the calculation for the same merger 
tree with altered numerical parameters as will be de- 
scribed below. The results of this exercise are shown in 
Fig. [5] in which we plot galaxies from this merger tree in 
the plane of total stellar mass and stellar mass-weighted 
scale length. The base model is shown by filled black 
circles. Galaxies in the comparison models which have 
a correspondent in the base model are indicated by open 
circles connected by a line to the base model galaxy. 
Where a base model galaxy is missing from the com- 
parison model it is marked with a cross, while in cases 
where a comparison model galaxy is not present in the 
base model it is indicated by a star. The parameters ad- 
justed in the three comparison models are: 

[odeToleranceRelative] Base model value: 0.01; 
Comparison model value: 0.001; Color: red. This 
parameter specifies the accuracy requested from 
the ODE solveij^] — it is required to keep fractional 
errors in evolved quantities below this value during 
any given time step. 

[timestepHostRelative] Base model value: 0.1; 
Comparison model value: 0.01; Color: green. 
This parameter limits the time interval by which 
satellite nodes may be evolved beyond the time at 
which their host halo is currently located. Specifi- 
cally, no satellite node is allowed to evolve beyond 
[timestepHostRelative] times the cosmologi- 
cal expansion timescales, H~ x (f), at the time of the 
host halo, before it must wait for the host node to 
catch up. 

[timestepSimpleRelative] Base model value: 0.1; 
Comparison model value: 0.01; Color: blue. This 
parameter, specified as a fraction of the cosmo- 
logical expansion timescale, H~ l (t), at the time of 
the node, limits the time step over which any node 
may be evolved before evolution is paused and any 
necessary post-processing (such as transferring gas 
driven out of a satellite galaxy) is performed. 

The results of this convergence study, as shown in 
Fig. [5] lead to several insights. First, the properties 
of half of the galaxies in this merger tree are almost 



8 We typically use a relatively large tolerance. Since the approxi- 
mations made by semi-analytic models are not expected to be valid to 
high precision a higher tolerance is not generally warranted. Higher 
tolerances can be adopted if necessary of course. 



23 



Table 1: Values of parameters used in the example model. Parameters selecting between different implementations of physical processes or 
components are only listed where more than one non-null implementation currently exists within Galacticus. 



Parameter 



Value 



Reference 



[H_0] 

[Qmega_0] 

[Omega_DE] 

[Qmega_b] 

[T_CMB] 

[accretionDisksMethod] 

[adaf Adiabaticlndex] 

[adaf EnergyOption] 

[adaf RadiativeEf f iciency] 

[adaf ViscosityOption] 

[adiabaticContractionGnedinA] 

[adiabaticContractionGnedinOmega] 

[barlnstabilityMethod] 

[blackHoleSeedMass] 

[blackHoleWindEff iciency] 

[bondiHoyleAccretionEnhancementHotHalo] 

[bondiHoyleAccretionEnhancementSpheroid] 

[bondiHoyleAccretionTemperatureSpheroid] 

[coolingFunctionMethod] 

[coolingTimeAvailableAgeFactor] 

[coolingTimeSimpleDegreesOf Freedom] 

[darkMatterProf ileMethod] 

[darkMatterProf ileMinimumConcentration] 

[diskOutf lowExponent] 

[diskOutf lowVelocity] 

[ef f ectiveNumberNeutrinos] 

[galacticStructureRadiusSolverMethod] 

[haloMassFunctionMethod] 

[haloSpinDistributionMethod] 

[hotHaloOutf lowReturnRate] 

[imf SalpeterRecycledlnstantaneous] 

[imf SalpeterYieldlnstantaneous] 

[imf Select ionFixed] 

[isothermalCoreRadiusOverVirialRadius] 



70.2 km/s 

0.2725 

0.7275 

0.0455 

2.72548 K 

ADAF 

1.444 

pure ADAF 

0.01 

fit 

0.8 

0.77 

ELN 

100 

0.001 

1 

1 

100 

atomic CIE Cloudy 


3 

NFW 

4 

2 

200 km/s 
4.34 

adiabatic 

Tinkei-2008 

Bett2007 

1.26 

0.39 

0.02 

Salpeter 

0.1 



54.2 



5 43 



5 43 



5 43 



§4.2 



5 43 



5 43 



5 43 



5 43 



5 43 



5 4.8 



5 43 



{VI 



533 



5333 



{3.1.2 



5333 



d Komatsuet al. || 2010[ ) 
(komatsu eTaTl feUTU I 



(KomatsuetaT 2010) 



( [Komatsu eTaT] [2TJTU I 
( [Komatsu et al.||201 o| 



5333 



5 433 



5 433 



5433 



5433 



5333 



5433 



5433 



5 433 



5 4.8 



5 43 



5 433 



5 333 



5 4332 



5 4.12.2 



54323 



5430 
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Table 1: (cont.) Values of parameters used in the example model. Parameters selecting between different implementations of physical processes or 
components are only listed where more than one non-null implementation currently exists within Galacticus. 



Parameter 



Value 



Reference 



[maj orMergerMassRatio] 

[mergerRemnantSizeOrbitalEnergy] 

[mergerTreeBuildCole2000AccretionLimit] 

[mergerTreeBuildCole2000MassResolution] 

[mergerTreeBuildCole2000MergeProbability] 

[mergerTreeConstructMethod] 

[minorMergerGasMovesTo] 

[modif iedPressSchechterFirstOrder Accuracy] 

[modif iedPressSchechterGO] 

[modif iedPressSchechterGammal] 

[modif iedPressSchechterGamma2] 

[powerSpectrumlndex] 

[power SpectrumRef erenceWavenumber] 

[power SpectrumRunning] 

[randomSpinResetMassFactor] 

[reionizationSuppressionRedshif t] 

[reionizationSuppressionVelocity] 

[satelliteMergingMethod] 

[sigma_8] 

[spheroidEnergeticOutf lowMassRate] 

[spheroidOutf lowExponent] 

[spheroidOutf lowVelocity] 

[spinDistributionBett2007Alpha] 

[spinDistributionBett2007LambdaO] 

[stabilityThresholdGaseous] 

[stabilityThresholdStellar] 

[starFormationDiskEf f iciency] 

[starFormationDiskMinimumTimescale] 

[starFormationDiskVelocityExponent] 

[starFormationSpheroidEff iciency] 

[starFormationSpheroidMinimumTimescale] 

[starveSatellites] 

[stellarPopulationPropertiesMethod] 
[summedNeutrinoMasses] 
[transf erFunctionMethod] 
[virialDensityContrastMethod] 



0.1 
1 

0.1 

5 x 10 9 M G 

0.1 

build 

spheroid 

0.1 

0.57 

0.38 

-0.01 

0.961 

1 Mpc- 1 



2 

9 

30km/s 
Jiang2008 
0.807 
1 

2 

50km/s 

2.509 

0.04326 

0.9 

1.1 

0.01 

0.001 Gyr 

-1.5 

0.1 

0.001 Gyr 
true 

instantaneous 


Eisenstein + Hu 
spherical top hat 
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Figure 5: Demonstration of convergence in galaxy stellar mass and 
stellar mass-weighted scale length. A single merger tree, with a final 
halo mass of W l2 M Q was run. Each point in the figure corresponds to 
a galaxy existing within that halo at z = 0. Black points correspond to 
our standard set of numerical parameters as defined in the text. Red 
points show a run with increased tolerance in the ODE solver. Green 
and blue points show runs with increased tolerance in time stepping 
relative to a node's host and to itself respectively. Where the same 
galaxy exists in the original model, the galaxy is represented by an 
open circle connected to the original model galaxy by a line. If a 
galaxy no longer exists in a modified model it is shown by a cross, 
while if a galaxy exists in the modified model that was not present in 
the original model it is shown by a star. 



perfectly stable to the specific changes in the numeri- 
cal parameters that we have explored, indicating a high 
level of convergence. The other three galaxies show sig- 
nificant changes in their properties however. We be- 
gin by considering the most massive galaxy. When 
the ODE solver tolerance is improved (red symbols), 
this central galaxy gains slightly in mass and becomes 
significantly larger. By tracing the merging history of 
this galaxy in Galacticus, we can identify the cause of 
this as a merger with the lowest mass galaxy that was 
present in the base model. In the base model calcu- 
lation, this lowest mass galaxy (located at a mass of 
around 3 x 10 6 M o in Pig. [5j» did not merge with the cen- 
tral. However, the improved ODE tolerance resulted in 
a merging event being triggered which then significantly 
altered the further evolution of the central galaxy. Con- 
versely, when the [timestepSimpleRelative] pa- 
rameter is reduced (blue symbols) we find a new galaxy 
appearing (indicated by the blue star). Here the oppo- 
site has happened — a galaxy which did merge in the 
base model failed to merge in the model with better 
timestep tolerance. This galaxy originally merged into 
the most massive galaxy in the tree. Consequently, that 
galaxy is now somewhat less massive and significantly 
smaller that in the base model. Finally, in this same 
model, the galaxy originally at 4.5 x 1O 8 M is increased 
in mass and slightly smaller. Once again examining the 
history of this galaxy in Galacticus, we find that this 
change occurs because the galaxy in question has a disk 
which sits close to the stability threshold for bar forma- 
tion (see §4.7| ). Under our current implementation of 
this process, the timescale for mass transfer from disk 
to spheroid due to the bar instability is finite at the sta- 
bility threshold, but infinite below that threshold. As 
such, properties of galaxies close to this threshold will 
be sensitive to small changes in numerical parameters in 
Galacticus. Fundamentally, this represents a limitation 
of the physical model adopted for the bar instability — 
clearly a more physically motivated model is desirable 
( Athanassoula 2008 ). While these effects make for sub- 
stantial changes in the properties of the galaxy popula- 
tion of this particular merger tree, we find much less 
change when considering statistical samples of galax- 
ies. Nevertheless, it is important to test for convergence 
in the quantities of interest for any particular study per- 
formed with Galacticus. 

Figure[6]shows another example of convergence. We 
take the same merger tree and re-simulate it while in- 
creasing the resolution of the tree branche^] Each line 



Since our trees are generated by a Monte Carlo method simply 
changing the resolution and rebuilding the tree would result in a very 
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in Fig.[6]represents a factor 5 decrease in the minimum 
mass resolved in merger tree branches. Solid lines show 
the cumulative stellar mass function of galaxies, while 
dashed lines show the cumulative mass function of dark 
matter (sub)halos scaled by a factor Ob/^o- Clearly the 
(sub)halo mass function is well converged. The galaxy 
mass function, however, is not. For example, the stel- 
lar mass of the most massive galaxy changes by a fac- 
tor of around 3 from the highest to lowest resolution 
runs. Of course, the results shown span a range of 625 
in merger tree mass resolution. This resolution depen- 
dence is expected — as lower mass halos are resolved 
more gas gets locked up in low mass galaxies leaving 
less to form the massive galaxy. Convergence should 
only be reached once a physical suppression scale is 
reached. In the case of cold dark matter this will be 
a baryonic scale (since any cut-off in the CDM power 
spectrum is expected to occur at very low masses) asso- 
ciated with a truncation in the cooling function or with 
the temperature of the IGM and photoheating by an ion- 
izing background. Even at the highest resolution shown, 
full convergence across the whole mass function is not 
reached. However, for resolutions of 2 x 1O 8 M and bet- 
ter the mass of the most massive (central) galaxy is con- 
verged to better than 50% which may be sufficient for 
many applications. Convergence in tree mass resolution 
will depend on the details of the implemented baryonic 
physics, and also on what galaxy samples/properties are 
being studied. This clearly demonstrates that a conver- 
gence test of the type carried out here should always be 
performed to ensure that results are not affected by the 
limited resolution of merger trees. 

5.1.2. Timing 

Galacticus was designed with simplicity, modular- 
ity and flexibility as key design principles. Execution 
speed was not a high priority but, nevertheless, we have 
given significant attention to optimizing key areas of the 
code. To explore the time taken to run typical calcula- 
tions we have evolved sets of merger trees of different 
masses in our example model and computed the mean 
time taken to evolve each tree. The results are shown 
in Fig. [7] where we show results for the full model, 
the same model but with black hole physics switched 
off and the same model but with all baryonic physics 
switched off. 



different tree (a single different branching would lead to an entirely 
different branching path). To make a fair comparison between the dif- 
ferent resolutions, we therefore create a tree at the highest resolution 
considered and then prune away branches below the required resolu- 
tion threshold. 




M. [MJ 



Figure 6: Tests of convergence in the cumulative stellar mass function 
of galaxies formed in a merger tree with a final halo mass of 10 12 M o 
at z = 0. Solid lines show the stellar mass function, while dashed 
lines indicate the (sub)halo mass function scaled by a factor of CI^/CIq. 
The difference in slope between the two sets of mass functions is a 
consequence of the SNe-driven feedback included in this model. It 
can be seen that the (sub)halo mass function is well converged, while 
the galaxy stellar mass function is not. Note that black holes were not 
included in calculations used for this particular convergence study. 
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Figure 7: Execution time per merger tree. (Run on a 2.7GHz AMD 
Opteron processor and compiled using standard optimizations.) The 
red region indicates the average time taken to evolve a tree in the ex- 
ample model to z = 0. The green region indicates the time taken to 
run the same tree with black hole physics switched off. Finally, the 
blue region shows the time taken to evolve each tree when only dark 
sector physics are included (i.e. all baryonic physics is switched off). 



Fig.|7]shows that the time taken to evolve a tree scales 
close to, but somewhat faster than, linearly with halo 
mass. Since we use a fixed mass resolution in these 
merger trees, the number of nodes into which a tree is 
resolved will scale approximately with the mass of the 
halo, so this result is to be expected. While a Milky 
Way-mass halo can be evolved in around 2 seconds, a 
high mass cluster takes around 1.25 hours. This is not 
as fast as some other semi-analytic models, but, as we 
stated above, speed is not our primary concern, with 
accuracy, simplicity and modularity taken precedence. 
Excluding black holes from the calculation significantly 
reduces run time by a factor of 2 to 2.5. This is because 
black hole properties (particularly spin) can evolve on 
very short timescales making evolution of the ODE sys- 
tem computationally expensive. This is worsened by 
the fact that many black holes are close to the equilib- 
rium spin at which spin-up by accretion and spin down 
from powering jets cancel, forcing small timesteps to 
maintain this balance. Finally, switching off all bary- 
onic physics (which leaves tree building, and evolution 
of the dark matter halos including dynamical friction, 
and which can be useful for exploring the properties of 
dark matter models) reduces run time by a further factor 



Figure 8: The z = bj-band luminosity function. Red points indicate 
the observed luminosity function from the 2dFGRS (Norb erg et al.| 
|2002| , while green points show results from our example Galacticus 
model after including the effects of dust-extinction using the model 
of |Ferrara et al.|fl999} . Error bars on the model points indicate the 
uncertainty due to the finite number of merger trees computed. 



of around 10. 

5.2. Galaxy Properties 

In this final subsection we compare several observa- 
tional datasets with their equivalents as predicted by the 
example Galacticus model. We reiterate that this ex- 
ample model is not intended to represent the best fit 
to observational data attainable with Galacticus — a full 
search of the Galacticus parameter space is a significant 
undertaking which will be explored in a future paper. 
Instead, the results in this subsection are intended to il- 
lustrate that types of quantities which can be computed 
and to show that they are, at least, in broad agreement 
with current observational constraints. We give a brief 
description and discussion of each result. Much of the 
underlying physics behind the model results has been 
previously discussed by other authors, as mentioned be- 
low. In a few cases, we defer further study to future pa- 
pers which will explore some of the key physical drivers 
in greater detail. 

Figures [8] and [9] show bj and K-band luminosity 
functions at z = derived from the 2dFGRS and 
2MASS galaxy surveys, with observational determina- 
tions shown by red points and results from the example 
Galacticus model show by green points. Dust extinc- 
tion is included in the model galaxy magnitudes using 
the results of |Ferrara et al.| (fl999). The match to the 
bright end cut off in the luminosity functions is quite 
good, although the K-band prediction undershoots the 
knee of the observed luminosity function. However, the 
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Figure 9: The z = K-band luminosity function from the exam- 
ple model. Red points indicate data from the 2dFGRS+2MASS 
( Cole et al. 2001 1, while green points show results from our example 
Galacticus model after including the effects of dust-extinction using 
the model of Ferrara et al. 1 1999). Error bars on the model points 
indicate the uncertainty due to the finite number of merger trees com- 
puted. 



predicted faint end slope is significantly steeper than 
those observed (particularly in the K-band). This is a 
well known problem in semi-analytic models (e.g. |Ben 



son et al.|20 03 ) and is closely linked to the inclusion of 
SNe-feedback in such models. Stronger feedback would 
help to reconcile this discrepancy. 

A key ingredient in producing the sharp cut off at the 
bright end of the galaxy luminosity function is the in- 
clusion of active galactic nuclei (AGN) feedback in the 



example Galacticus model ( Croton et al. |2006| |Bower| 
|et al.|[2006| . As such, it is crucial that predicted black 
hole properties be reasonable. Figure [10] shows the rela- 
tion between black hole mass and spheroid stellar mass 
from observations (red points) and from the example 
Galacticus model (green points with error bars; indicat- 
ing the mean and l-<x dispersion at each spheroid mass). 
The normalization and slope of the model relation is in 
good agreement with that which is observed. The model 
displays little scatter in black hole mass at fixed stellar 
mass, except for at the highest masses (partially due to 
the limited number of such high mass spheroids in our 
model sample) and around 1O 1O M spheroid mass. In 
the model, this increase in scatter below 10 10 M o oc- 
curs because the regulating mechanism of AGN feed- 
back (which controls the relation between black hole 
and spheroid mass) breaks down. 

Figure [TT] explores the gas content of galaxies in the 
example model, by comparing the model HI mass func- 
tion with that observed by |Zwaan et al.| ( f2005| l. While 
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Figure 10: The relation between supermassive black hole mass and 
galaxy bulge stellar mass. Data are taken from Haring and Rix (2004) 
and are shown by red points. Green points with error bars show the 
mean and 1-cr dispersion from our example Galacticus model. 
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Figure 1 1 : The z = HI gas mass function. The green points show the 
mass function from the example Galacticus model assuming a con- 
stant factor of 0.54 to convert from total gas mass to HI mass (Power| 
|et al.| [2010 l. Red points indicate data from Zwaan et al. 1 2 005} . 
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Figure 12: 01 g- 01 r differential color distribution (normalized to unit 
area under the surface) for galaxies at z = 0. 1 as a function g-band ab- 
solute magnitude. The color scale shows data from the SDSS jWein-| 
mann et aL 2006 1, while the black contours show the distribution from 
the example Galacticus model. Ten contour levels are shown, linearly 
spaced between the minimum and maximum of the color scale. 



the overall normalization is approximately correct the 
shape of the predicted function is incorrect and, in par- 
ticular, does not truncate sharply enough at high masses. 
Similar results were found and considered in greater de- 
tail by |Power et al^pOlO) . 

A key result from the SDSS has been a clear demon- 
stration of the dichotomy of the galaxy population, with 
most galaxies being part of a star forming "blue cloud" 
or a quiescent "red sequence". Figure [12] shows the 
SDSS color-magnitude diagram (colored shading) de- 
rived by Weinmann et al. ( 2006) > which clearly shows 
these two populations. The black contours show the 
color-magnitude diagram from the example Galacticus 
model. This clearly shows the same two populations of 
galaxies, with approximately the correct median colors 
and transition luminosity. In the example model, AGN 
feedback plays a key role in establishing this division, 
as has been found by other authors ([Croton et al.|2006 



|Bower et al.|2006| see also |Font et ai.|2008) . 



The structural and dynamical state of model galaxies 
is explored in Figure [13] in which we compare results 
from the example model to the i-band Tully-Fisher rela- 



tion measured by Pizagno et al. ( 2007 1. The model is in 
moderately good agreement with the data — the slope at 
lower luminosities is approximately correct, but flattens 
at higher luminosities such that the model velocities are 
too low. Additionally, the dispersion in model galaxy 
velocities at intermediate luminosities is too large. The 
Tully-Fisher relation encodes significant constraints on 
the process of galaxy formation. As such, we intend to 
explore this aspect of the Galacticus model in greater 
detail in a future paper. 



Pizagno et al. (2007) > 
Galacticus (mean+dispersion) 




Figure 13: The i-band Tully-Fisher relation from the SDSS i Pizagno 
|et al.| [20 07 1 is shown by red points, while results from the example 
Galacticus model are indicated by green points. Model points show 
the mean disk circular velocity at each magnitude with error bars in- 
dicating the l-cr dispersion in model galaxies. 




Figure 14: The star formation rate per unit comoving volume in the 
Universe as a function of redshift. Red points show observational 
estimates from a variety of sources as compiled by |Hopkins| {2004} 
while green points show the star formation rate inferred from gamma 
ray bursts by Kistler et al. C2009). The blue line shows the result from 
the example Galacticus model. 
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Figure 15: Distribution of disk scale lengths for galaxies at z = se- 
lected by face-on I-band absolute magnitude, -20 < M\ q < -19. Red 
circles show data from de Jong and Lacey 1 2000 1 with upper limits in- 
dicated by red triangles. Green points show results from our example 
Galacticus model. 



The volume averaged star formation rate density in 



the example model is shown in Fig. 14 and is compared 
with a compilation of observational data. Clearly the 
model under-predicts the star formation rate at all red- 
shifts above z ~ 0.25. This is a key failing of the ex- 
ample model, and therefore represents a key constraint 
for future parameter space studies with Galacticus. The 
origins of this underestimate of the star formation rate 
are not clear without further investigation, but could 
plausibly be due to over-zealous AGN feedback at in- 
termediate redshifts or simply a breakdown in the as- 
sumed scalings of star formation timescales (which are 
currently based on simple scaling relations rather than 
any physical model). 

Finally, we show in Fig.[15]the sizes of galactic disks 
for galaxies in a narrow range of face-on I-band mag- 
nitude, -20 < Mi o < -19, and compared to the ob- 
servational determination of |de Jong and Lacey (2000). 
The peak of the distribution in the example model is in 
approximately the correct location, but the model pre- 
dicts somewhat too many very large disks compared to 
the data (although, given the large errors on the data in 
this region, the difference between model and data is not 
very significant). Sizes of galaxies are a key property 
since they directly affect galactic dynamical timescales 
and, therefore, star formation rates. 

6. Summary 



We have described a new, free and open 
source semi-analytic model of galaxy forma- 



tion, Galacticus, which can be downloaded from 
http : //sites . google . com/site/galacticusmodel 
A full manual documenting both the physics and tech- 
nical implementation of Galacticus can also be found 
on the same websit^] In particular in this article, 
we have focused on the technical and practical im- 
plementation of Galacticus and have detailed the 
numerous physical ingredients which make up the 
model. Additionally, we have explored basic properties 
of a specific example model, demonstrating numerical 
convergence and illustrating the type of quantities 
which Galacticus can compute and compare against 
observational data. Galacticus contains the majority of 
the physics incorporated into other major semi-analytic 
models (see |Benson| ( |2010| l for a review). Missing 
physics, notably tidal and ram pressure mass loss from 
satellite galaxies (|Benson et al.| |2002| |Lanzoni et al 



2005 Font et al. |2008t |Henriques and Thomas] |2010 



and the ability to use information on subhalos from 
N-body merger treesj^] is already being incorporated 
into the development version of Galacticus. 

The Galacticus model can be employed for a wide 
variety of calculations related to galaxy formation stud- 
ies, ranging from dark matter phenomenology, through 
observationally testable statistical properties of galaxy 
samples to the construction of mock galaxy samples. 
Its flexibility makes it easily adaptable to new problems 
and accepting of new physical prescriptions. 

Galacticus requires only free and open source li- 
braries and compilers to run and so is portable and 
practical to use. In designing Galacticus we placed 
an emphasis on simplicity, flexibility and extensibility 
to make it simple to use and maintain and to facilitate 
future development as required in the rapidly chang- 
ing field of galaxy formation. We therefore hope that 
Galacticus will prove to be a valuable tool for a signif- 
icant time, and we are already endeavoring to expand 
and improve upon this initial version. 
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10 This article specifically describes the current development ver- 
sion of Galacticus, vO.9.0. The current stable version, v0.0. 1, con- 
tains almost all of the same features. 

"Galacticus can utilize simple merger trees from N-body 
simulations — software is provided, for example, to convert trees from 
the Millennium Simulation into Galacticus's format — but not those 
with subhalo information. 
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